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Abstract. Low Reynolds number direct simulations of large populations of 
hydrodynamically interacting swimming particles confined between planar walls are 
performed. The results of simulations are compared with a theory that describes 
dilute suspensions of swimmers. The theory yields scalings with concentration for 
diffusivities and velocity fluctuations as well as a prediction of the fluid velocity 
spatial autocorrelation function. Even for uncorrelated swimmers, the theory predicts 
anticorrelations between nearby fluid elements that correspond to vortex-like swirling 
motions in the fluid with length scale set by the size of a swimmer and the slit height. 
Very similar results arise from the full simulations indicating either that correlated 
motion of the swimmers is not significant at the concentrations considered or that 
the fluid phase autocorrelation is not a sensitive measure of the correlated motion. 
This result is in stark contrast with results from unconfined systems, for which the 
fluid autocorrelation captures large-scale collective fluid structures. The additional 
length scale (screening length) introduced by the confinement seems to prevent these 
large-scale structures from forming. 



Submitted to: J. Phys.: Condens. Matter 



Dynamics of confined suspensions of swimming particles 



2 



1. Introduction 

Experimental observations of suspensions of swimming microorganisms illustrate a 
number of fascinating phenomena that are still poorly understood. Correlations of 
the swimmers result in jets and swirling motions on scales larger than that of a single 
organism [K [2j. The organisms also form local nematic ordering though they have no 
global nematic behavior [3j. The collective behavior leads to swimmer velocities larger 
than that of an isolated organism [4j and enhanced transport in the fluid [5j. 

There are many open questions regarding these important observations both in 
terms of what leads to the phenomena and their biological signiflcance. Some models 
consider "local" interactions between nearby organisms within a flnite range, and how 
those can lead to large-scale collective behavior. These interactions may take the form 
of ad-hoc rules [6j or direct steric interactions between the swimmers Other models 
consider long-ranged hydrodynamic interactions between swimmers, which decay as 
in an unbounded domain (HJ |9l [IHl El [12] . The relative importance of these phenomena 
is still not fully understood. 

Many experiments observing these phenomena have been performed in droplets or 
thin fllms, but the influence of conflnement on the observations is not clear. Conflnement 
sterically hinders the organisms and also aflFects the hydrodynamic interactions, both 
between the organisms and between the organism and the boundary. The boundaries 
may also play a role in transport of nutrients, which may, in turn, affect motion of 
the organisms. For example, oxygen levels may be different near surfaces, altering the 
motion of the organisms. For droplets or fllms, the rigidity of the surface from secreted 
molecules may affect the fluid boundary condition and dynamics of the organisms [H [3] . 

Most experiments mentioned above have concentrated on two types of bacteria, 
E. coli and B. suhtilis. Both organisms use flagella to propel themselves forward through 
the fluid. However, other types of organisms propel through a fluid by other mechanisms, 
such as pulling themselves forward from the front ^13j . What role the propulsion method 
plays on the collective behavior has not been clarifled, although computational and 
theoretical evidence suggests that organisms pushed from behind display more complex 
collective behavior than those pulled from the front [Ml [12]. In general, what biological 
signiflcance the above phenomena play in the function of the microorganisms remains 
unclear. 

In this article we focus on the role that conflnement plays in the collective behavior. 
The only prior publication using computational models to investigate how hydrodynamic 
interactions affect the collective behavior of swimming microorganisms in conflned 
environments is by Hernandez et al. [8j. They showed that hydrodynamic interactions 
were sufficient to produce many qualitative phenomena seen in experiments, including 
increased transport and swirls in the fluid. We focus on conflned systems not only 
because many experiments have been done in conflned environments, but also to draw 
comparisons with unconflned dynamics. Simulations in three-dimensional periodic 
domains have shown changes in dynamics with the size of the domain [12j that are 
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thought to arise because of long-range structures that fiU the entire domain. Some 
experiments have attempted to measure unconfined dynamics by examining the response 
far away from the confining waUs p!5l [T6l \T7\ . However, because of the long-ranged 
nature of hydrodynamics, the point at which confined systems with very large gaps 
reduces to an unconfined system is not obvious. 

We examine the behavior at diflFerent levels of confinement and at different 
concentrations ranging from the dilute limit into the semidilute regime. However the 
concentration is still low enough that we do not expect the steric interactions between the 
organisms to be the dominant interaction and lead to nematic-like structures. Instead 
we expect the long-ranged hydrodynamic interactions to play an important role. The 
confinement boundaries alter these hydrodynamic interactions. The walls also induce a 
non-uniform concentration profile within the domain, with the organisms concentrating 
at the walls. Finally, the walls introduce a new length scale which screens hydrodynamics 
and alters flow structures over larger length scales. We will show in this article how the 
dynamics depend on conflnement and how they compare with theoretical predictions. 

2. Swimmer Model 

Consider a suspension of neutrally buoyant rodlike swimmers confined between two 
planar walls. Directions Xi and X2 are periodic of side length L, and the walls are 
separated in the X3-direction by a distance 2H. Each swimmer has a characteristic 
length ^, a characteristic width and in isolation would move in a straight line with 
a speed Vi^. It is assumed that the Reynolds number. Re = v^^H/v <C 1, where v is 
the fiuid kinematic viscosity, in which case the fiuid motion is governed by the Stokes 
equation. To allow treatment of large populations (> 10^ swimmers) over long times, a 
very simple model of each swimmer is adopted, following our previous work jHl ri2] . 

Each self-propelled particle is modeled as two beads connected by a stiflF spring 
with equilibrium length I as shown in figure 1. The unit vector pointing along the 
swimmer axis from bead 1 to bead 2 is denoted n. Propulsion is provided by a "phantom 
fiagellum" that we do not treat explicitly, but only through its effect on the swimmer 
body and the fiuid. The bead which is connected to the fiagellum, denoted bead 1, 
feels a force exerted on it. However, the fiagellum also exerts a force — on the 
fluid. This force on the fluid occurs at the position of bead 1. With this model we 
can consider "pushers" or "pullers" depending on whether is parallel or anti-parallel 
to n, respectively. A pusher sends fluid away from it fore and aft, with fluid moving 
toward its "waist", and vice versa for a puller. Whether a real swimmer is a pusher or a 
puller depends on the speciflc mechanism of locomotion (cf. p^)- a cell whose flagella 
propel it forward predominantly from behind would be a pusher [HI [12]. We will focus 
primarily on pushers in this article, though some results with pullers are presented. 

To measure the concentration of the swimmers we deflne an effective volume fraction 
(j)Q = TiNi^ /{12L'^H) - this would be the true volume fraction if the swimmers were 
spheres of diameter i. In the present geometry, swimmers in dilute systems form layers 
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Figure 1. Illustration of a pushing organism, our swimmer model, and the fluid 
disturbance they cause in an unbounded domain. Our model shows the hydrodynamic 
radius of the beads and the ellipsoidal excluded volume. The double-arrow signifies the 
phantom fiagellum force acting on the bead and an opposite force acting on the fiuid. 
Both forces act at the center of the first bead. The blue lines represent streamlines of 
the axisymmetric fiuid disturbance. A puller would produce the same streamlines, but 
with the direction of the velocities reversed. 



near the two walls. Therefore, we find it convenient to also define an effective area 
fraction if all swimmers resided in the layers. It is defined as ip^ = CnNP /{4:L'^)^ 
where the constant (7 = 1/2 if the gap is large enough to allow for a layer at each wall 
{2H > 2w)^ or (7 = 1 for extreme confinements, e.g. monolayers. 

The force in the "stiff" connector between the beads follows a finitely extensible 
non-linear elastic (FENE) spring model with a non-zero equilibrium length equal to the 
swimmer size £ [18j: 

where h is the spring constant, is the connector vector from bead 1 to bead 2 on 
swimmer z/, and £^ is the swimmer maximum size. With this spring law, the swimmer 
will shrink slightly for a pushing fiagellum, or expand for a pulling one. The values of h 
and are chosen so that the spring approximates a rigid constraint while still allowing 
timesteps that are not prohibitively small. For the results presented here, h = f /{0.1£) 
and = 1.15^. The swimmers also interact through an excluded volume potential, 
which is taken as the repulsive portion of the Gay-Berne potential [19] ; this potential is 
widely used in molecular simulations to model steric repulsions between rodlike objects. 
The size and aspect ratio of the excluded volume potential to related to the size of the 
swimmer £ and the hydrodynamic radius of a bead a as illustrated in figure [2] The 
width of the excluded volume is taken slS w = 2a and the length as ^ + 2a. This means 
that the excluded volume restricts bead positions such that the hydrodynamic radii do 
not overlap. It also gives an aspect ratio of 7 = 1 + £/{2a). The model presented here 
has ^ = 3a, and therefore an aspect ratio of 7 = 2.5. 

The motion of the swimmers is determined by the force balance (neglecting inertia 
because of the small size of a microorganism) for each bead (fc = 1, 2) of a swimmer z/. 
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as follows 

t^ + t% + f^fc + 4if: = for k = l,2, (2) 

where 5ij is the Kronecker delta, i^j^ is the hydro dynamic drag force, i^j^ is the connector 
(spring) force and f^^ are the bead-bead, bead-swimmer and bead- wall excluded volume 
forces. Notice that the force balance for bead 1 differs from the balance on bead 2 by 
the presence of the flagellum force, f^. 

In our simulations, each bead will be treated as a point particle. In this situation, 
the hydrodynamic drag force on a bead k of swimmer u is given by a generalization of 
Stokes' law [20j: 

fjfc = -C (v^,fc - u^,fc) , (3) 
where ( = Gnrja is the Stokes drag coefficient on a bead with hydrodynamic radius a 
in a fluid with viscosity ry, v^^/^ = ^iy,k is the velocity of the bead, where x^^^^^ is the 
position (cartesian coordinates) of the bead, and Uj^^k is the fluid velocity at the bead 
position. There are two contributions to this fluid velocity: the flrst is the motion driven 
by the forces exerted by the other beads (and flagella) in the system, and the second is 
the correction to the velocity experienced by the bead due to the presence of conflning 
walls. This correction leads for example to the decrease in bead mobility found in a 
conflning geometry. Both of these contributions are determined simultaneously by the 
methodology used here. Equations ^ and ^ give an evolution equation for the bead 
positions as follows, 

^ = u.. + ^ (t% + C,. + . (4) 

For i?e = 0, the fluid velocity at a point x due to a collection of point-forces 
is calculated by summing over all forces times the Green's function G for the Stokes 
equation for the geometry and boundary conditions of interest. Considering the beads 
of each of the swimmers as point-forces and including the disturbance due to the 
phantom flagellum, we write the fluid velocity at a point x as 

2 

^ W = E E x,,0 • (-fi, - Sail) , (5) 

1=1 

where the flrst term represents the direct hydrodynamic forces exerted by the beads 
on the fluid and the second term represents the disturbance caused by the phantom 
flagella. These two contributions can be combined using the force balance on each bead 
to become 

N 2 

^(^) = E E G(x, x^,o • (f;,, + . (6) 

Similarly, the fluid velocity at the position of a bead k of swimmer u is calculated by 
excluding the singular part of the fluid velocity generated by the bead: 

N 2 

12=1 1=1 
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Here Goo is the Oseen-Burgers or Stokeslet tensor, the free-space point-force Green's 
function for Stokes equation. It is given exphcitly as: 

with r = |x| and 6 the identity tensor. In an unbounded domain, G and Goo are 
identical so the u = /a^ k = I terms in this expression are zero. In a bounded domain, 
however, G contains a finite correction due to the change in the fiow induced by the 
boundaries; subtracting off Goo when u = /a^ k = I reveals this finite correction in 
the overall expression. Substitution of Eq. ([7| into the evolution equation for the bead 
positions, Eq. ([i]), results in 

^ 1=1 

Here M(^^/e)(/i,o is a 3 x 3 tensor which constitutes a block of the (3 x 2A^) x (3 x 2A^) 
mobility tensor, M: j2T] 

^{u,i){ii,k) = Siy^5ki- + (G(x^,fc, x^,^) - 5^^4^Goo(x^,/c - x^,/))- (10) 

Equation ^ can be written in a compact form by introducing 3 x 2N dimensional 
vectors containing the coordinates and forces of all beads as follows: 

^R= ^F^ + M-F^, (11) 
dt C 

where R contains the coordinates of all beads, F^ is a vector with the fiagellum 
force, whose components are non-zero for bead 1 of each swimmer, and F^ is a vector 
containing the total non-hydrodynamic and non-fiagellar forces on each bead. 

The net force exerted by an isolated swimmer on the fiuid is zero - to leading order 
in the far field a neutrally buoyant swimmer is a force dipole. The present approach 
captures this universal far-field behavior while neglecting the near-field corrections to 
the hydrodynamic interactions between swimmers, which are dependent on the details 
of the organism. The validity of this approximation is supported by recent simulation 
results fTT] . Finally, we note that the limited set of results we have obtained with multi- 
bead rod swimmers is qualitatively consistent with those for the two-bead swimmers. 

The fiuid velocity M • F^ is calculated using the General Geometry Ewald-like 
Method (GGEM) introduced by Hernandez-Ortiz et al [22]. A brief description of the 
GGEM starts with considering the Stokes system of equations for a fiow driven by a 
distribution of 2N point forces, 

- Vp(x) + ryV Vx) = -p(x) 

V.u(x) = ' ^ ' 

where rj is the fiuid viscosity and the force density is p(x) = Ylt=i fi^^(x — x^). Here ij^ 
is the force exerted on the fiuid at point x^. The solution of ( [T2| ), can be written in the 
form of ([g]) and combined into the M • F product. If computed explicitly, this product 
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is a matrix-vector operation that requires 0{N^) calculations. GGEM determines the 
product implicitly for any geometry (with appropriate boundary conditions) without 
performing the matrix- vector manipulations. It starts with the restatement of the force- 



density expression in (12), p(x) = Pi(x) + Pg(x) using a smoothing function similar 
to conventional particle-mesh methods. By linearity of the Stokes equation, the fluid 
velocity is written as a sum of two parts, the solution with each force-density separately. 
The "local density" 

2N 

drives a local velocity, Ui(x), which is calculated assuming an unbounded domain: 
Ui(x) = Gi(x — x^) • f^, where Gi(x) is composed of the free-space Green's function 
Goo minus a smoothed Stokeslet obtained from the solution of Stokes equations with 
the forcing term modifled by the smoothing function g. For the Stokes equations we 
found that a modifled Gaussian smoothing function deflned by 

g{r) = (aV7r^/')e(-"'"'H5/2 - c^V') (14) 
yields a simple expression for Gi(x): 

Wc xx^ 2a ^(_^2,2) 



Because Gi(x) decays exponentially on the length scale in practice the local velocity 
can be computed, as in normal Ewald methods, by only considering near-neighbors to 
each particle u [23j. 

For the present work, the point-particle approximation is not desired; in particular, 
as the concentration of particles increases, the probability that particles will overlap, 
having unphysical velocities, increases. To avoid this problem, the bead hydrodynamic 
radius, a, can be used to deflne a new smoothed-force density which will give a non- 
singular velocity. This is achieved by replacing the Stokeslet by a regularized Stokeslet, 
using the same modifled Gaussian with a replaced by with ^ ^ a~^, yielding 

erf(^r) eii^ar) 



^ ^ ^ 87r?7 V r2 y 



where the superscript R stands for regularized force density. For = 3a /^/n^ the 
maximum fluid velocity is equal to that of a particle with radius a and the pair mobility 
remains positive-deflnite [221 [24] . 

The global velocity, Ug(x), is due to the force distribution Pg(x), which is given by 

2N 

^g(^) = ^^iugi^-^u)' (17) 



iy=l 
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For a general domain, we find the solution to Stokes' equation numerically, requiring that 
Ui + Ug satisfy appropriate boundary conditions. At a no-slip boundary we would require 
Ug(x) = — Ui(x). For problems with periodic boundary conditions, Fourier techniques 
can be used to guarantee the periodicity of the global velocity Ug. The periodicity on 
the local velocity, Ui, is obtained using the minimum image convention. In the present 
case, the global contribution is calculated on a mesh with M = M1M2M3 mesh points, 
with Ml, M2, and M3 the number of mesh points in the Xi, X2, and Xs directions. A 
Fast Fourier Transform (FFT) method is used in the two periodic directions {xi and X2) 
and a second order finite difference method (FDM) scheme for the confined direction 
(xs). The FFT is implemented using FFTW ^ El, while the FDM is solved with a 
regular LU decomposition routine [27] . 

For GGEM, the results should be independent of a but the computational cost of 
the local and global calculations depend on a. Therefore, we choose the optimal a which 
minimizes the total computational cost. In the global calculation, to have an accurate 
solution, the mesh size must be smaller than the scale of the smoothing function, which 
is a~^. Therefore, Mi^2,3 ^ oi. The cost of the each LU decomposition scales as M| 
and there are M1M2 number of them, giving a total global cost that scales as . Note 
that the cost of each FFT scales as MiM2ln(MiM2) and there are M3 of them giving 
a scaling of a^ln(a). This is small compared to the cost of the LU decomposition. In 
the local calculation, the contribution of all pairs must be calculated that lie within a 
neighbor list determined by the decay of the local Green's function. The local Green's 
function decays over a distance so the number of neighbors for each particle scales 
as Na~^. The calculation must be performed over all pairs, which is the number of 
particles times the number of neighbors per particle, resulting in a local calculation 
cost scaling of N^a~^ . Minimizing the total (local and global) computational cost with 
respect to a gives an optimal a that scales as a^^t ^ A^^/^ and a total cost that scales 
as (9(A^/^). If we had chosen a diflFerent, linear, method for the FDM calculation, the 
global cost would have scaled as a^, leading to an optimal value of aopt ^ A^/^ and a 
total computational cost that scaled as 0{N). 

For our slit geometry, we used for the periodic directions a number of mesh points 
^1,2 = V^o^L while M3 = 4:^/2aH for the confined direction. The previous analysis 
determined how the optimal value of a changes with system size A. However, to 
determine the value of the prefactor, and thus the value of a used in the simulations, 
the cost of the global and local contributions must be determined on the computers 
used for the computation. For our machines, determining the computational cost led to 
using ttopt = 0.042A1/1 

3. Theory for Dilute Suspensions 

As a starting point for understanding the dynamics of confined suspensions of swimmers, 
we present here a scaling theory valid in the dilute limit, in which the swimmers 
act almost independently. We also show the fiuid correlations for a suspension of 
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uncorrelated swimmers. These properties are an important baseline in order to identify 
whether or not an observed feature of a swimming suspension arises from coUective 
dynamics. 

3.1. Diffusivity and velocity scaling 

In the dilute limit, swimmers are mostly to be found in layers near each wall [8]. The 
fundamental reason for this is simple - a swimmer that is not near a wall is in all 
likelyhood oriented toward one wall or the other and will eventually collide with it. Once 
it collides with a wall it remains there until a fluctuation causes it to leave (in which 
case it will eventually hit the other wall) . We shall see that swimmer- wall hydrodynamic 
interactions have a quantitative effect on the layering, but not a qualitative one. Recent 
experiments with bacteria also show a very high concentration at the walls [28] . Because 
the formation of layers is a generic phenomenon, the present theory focuses on the 
dynamics of swimmers in the layers. 

In the absence of other swimmers, a swimmer will continue in a straight line within 
the layer. Collisions or hydrodynamic interactions with other swimmers cause the 
direction of swimming to change, leading to diffusive motion at long time scales. From 
a scaling perspective, consider a random walk in which the particle swims at constant 
speed along a trajectory for a mean duration before changing directions. The 
motion at long times is diffusive with a diffusivity that scales as ^ v^r^ ^ vj^^ where 
= ^s^s is the mean free path. In a dilute system, the velocity of the swimmer is 
essentially the isolated swimmer value, i.e. Vg ~ Vig. (Changes in swimming speed due 
to conflnement do not affect the scaling predictions.) Deviations from this at larger 
concentrations will be discussed later. The mean free path scales as /g ^ ^'^{^i^e)~^^ 
where a is a two-dimensional cross section (with units of length) for the redirections. 
Therefore, the swimmer diffusivity, Dg, will scale as 



in dilute systems. Analysis for an unconflned domain [12] yields a similar scaling but 
with a three-dimensional cross section ctsd and the area fraction replaced by a volume 
fraction, so A,3d ^ (crsD^e)"^- 

We now turn to the fluid flow generated by the motion of the swimmers, considering 
the behavior of passive, non-Brownian tracers that follow the local fluid velocity. The 
tracers undergo ballistic motion at short lag times and diffusive motion at large lag 
times. The velocity of a tracer, Vt, is the fluid velocity at the location of the tracer Xt, 
and is the sum of the disturbance velocities due to each swimmer: 



where is the disturbance at the position of the tracer due to swimmer u. We can 
calculate the mean-squared velocity of the tracers by squaring this sum and performing 



-1 



(18) 



N 




(19) 
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an ensemble average over all possible configurations of the swimmers while the position 
of the tracer is held fixed at Xt: 



N N 

(vfvO.. = K>., = EE«-<>x.' (20) 

where the subscript Xt denotes that the tracer is held fixed at Xt during the ensemble 
average over swimmers. The ensemble average without such a subscript means an 
average over all tracer positions has also been performed. In the dilute limit, we 
assume the swimmers to be distributed independently. This means that each swimmer 
is independently sampling a probability distribution for location within the channel, 
though that probability distribution need not be uniform. For example, we consider 
here that distribution to be peaked in a layer near each wall. The independence of the 
swimmers leads to 

K>.,-A^(W)^^-^e^^e. (21) 

The number of swimmers is proportional to both the volume fraction and the area 
fraction with a different proportionality factor that depends on the channel height. 

A similar scaling analysis can be used for the tracer diffusivity, which can be written 
using the Green-Kubo relation [29l [30] 



Dt = i, (vt(0)-vt(t))dt (22) 



Again, we can replace the tracer velocity by a sum over swimmer disturbances. Assuming 
independent swimmers in the dilute limit gives 

A ~ y <u^(0) ■ utit)) dt - (/>e ~ A. (23) 

Combining the velocity and diflFusivity scalings with the definition of a correlation time, 
Tt = Dt/{vl) gives a prediction that is independent of volume fraction in the dilute 
limit. 

We can also use this behavior of the tracers to understand the correction to the 
swimmer behavior at finite concentration. The tracers represent a characteristic fiuid 
element in the system. In addition to being self-propelled, each swimmer is advected 
by the local fiuid disturbance (generated by the other swimmers). If the fiuid velocity 
does not change rapidly over the size of a swimmer, then we can write 

Vs ^ Vis + Vt(Xs). (24) 

Squaring this expression and neglecting the cross terms in the dilute limit gives 

<-s^> - 4 - {vl) . (25) 

Note that there exists a proportionality factor because the average over the fiuid in the 
region occupied by the swimmers is different than the average over the fiuid in the whole 
domain. The proportionality is replaced by an equality for an unconfined system, for 
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which these domains are the same [12j. Combining this result with the scahng of the 
tracer velocity {v'^) ^ (j)^ ^ ip^ gives that: 

{Vl) - Vi r^<P,r^ (26) 

and similarly that 

3.2. Velocity fields and spatial correlations 

An important measure of correlations in a suspension of swimming microorganisms 
is the spatial autocorrelation function of the fluid velocity. We will report results 
for correlations on planes of constant Xs- The separation of different fluid elements 
in this plane is denoted by xy. Therefore, the fluid correlation is (7f(x||,X3) = 
(u(s||, Xs) • u(s|| + X||, X3)), where the angle bracket represents an ensemble average over 
swimmers and also an average over the tracer position sy in the plane X3 = const. 

A number of experimental studies [21 |3j have reported the spatial correlation 
function for the swimmer velocities: (7s(x||,X3) = (vs(s||, X3) • Vs(s|| + xy, X3)). 
Assuming the validity of the approximation given by Eq. [24| we can write that 

Cs{^\\,Xs) = (vis(S||,X3) • Vis(S|| +X||,X3)) 

+ 2(vis(S||,X3) •U(S|| +X||,X3)) (28) 
+ Cf(x||,X3). 

In the case of independent swimmers whose orientations are uncorrelated with the fluid 
velocity, this expression reduces to 

Cs(x||,X3) =^fgXo(X||)+Cf(x||,X3), (29) 

where Xo(x||) = 1 if xy = and zero otherwise. This expression indicates the close 
relationship between Cg and Q. In our simulations, Q is substantially less susceptible 
to statistical noise than Q so that is the quantity we report. 

The correlation function Cg has been used in experiments and simulations to 
quantify the "swirls" seen in the fluid |21 |3j. In particular, the existence of negative 
correlations has been used as evidence for collective behavior and used to quantify the 
size of the swirls. However, we will see here that negative correlations can be present even 
in the absence of collective behavior. The presence of the walls changes the disturbance 
that an organism produces in the fluid, and thus the correlations in the fluid. 

The spatial fluid correlations are determined by the fluid disturbance caused by a 
single swimmer as well as correlations between the swimmers. We consider here the 
correlations present in the absence of correlations between the swimmers, that is for 
independent swimmers. The correlation for independent swimmers is determined solely 
by the disturbance caused by a single swimmer and the concentration of swimmers. 
Recall from flgure 1 the disturbance that a single swimmer produces in the fluid in the 
absence of walls, a pusher expels fluid out from the front and back, while sucking in fluid 
from the sides. However, the streamlines do not form closed loops, or swirls. Contrast 
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Figure 2. Streamlines of the fluid disturbance due to a force dipole representing a 
pusher: at the middle of the channel (top row) and close to the wall (bottom row). 
The fluid is shown in the X1-X2 plane (left column) and the Xi-Xs plane (right column). 

this with figure 2, which shows the streamlines of the disturbance produced by a single 
swimmer in the presence of waUs. The waUs induce swirls in the flow. Two new length 
scales, the separation of the walls and the separation of a swimmer from the walls, affect 
the fluid structures. 

The difference between these flow flelds leads to different fluid correlations. As we 
show below, the swirl in the fluid induced by the wall means that even independent 
swimmers, that have no collective behavior, will produce fluid correlations that are 
negative. To illustrate this point, it is useful to visualize snapshots of a typical velocity 
field generated by these independent swimmers, as shown in figure |3) It is important 
to note that these snapshots are simply a sum of disturbances due to independent 
swimmers, each producing a disturbance like that in figure 2. The snapshots look 
remarkably similar to the results of simulations and experiments for which the swirls 
were considered evidence of collective behavior. 

The autocorrelation function within a fiuid for which there are such swirls has a 
region of negative correlation, which is shown in figure |4} The correlation function has 
been calculated in two planes, the plane in which the swimmers form layers as well 
as the center plane of the slit. By changing the amount of confinement, we see that 
the length scale at which the correlation becomes negative in the center of the slit is 
the separation of the walls. However, the correlation near the wall, in the plane of the 
swimmers, is independent of the confinement. Instead, the length scale depends on the 
size of a swimmer and the separation of the swimmers from the wall. This is true for 
the cases shown because the second wall is far enough away that its effect on the fiuid 
correlation is small. 
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Figure 4. Fluid velocity autocorrelation function for independent swimmers at 
2H = 5i and 2H = 10£. The autocorrelation is normalized by the corresponding 
value at |x||| = 0. When xs = H - £/2 the curves for 2H = U and 2H = 10£ are 
indistinguishable. 



By examining theoretically the fluid correlations generated by uncorrelated 
swimmers, we have shown that collective behavior is not necessary to generate swirls 
and negative correlations in the fluid. This illustrates the importance of comparing 
properties calculated from a suspension of swimming organisms to that of independent 
organisms to gauge the importance of collective behavior in producing the observed 
response. 

4. Computational Results 

The theory just described makes distinct predictions regarding swimmer and tracer 
diffusivities, velocities, and fluid phase correlations. In this section we directly compute 
these quantities and others for suspensions over a range of concentrations and degrees 
of conflnement. For the remainder of the article, all lengths are non-dimensionalized by 
the swimmer size £ and all times are non-dimensionalized by i/vis- 

Let us start by flxing the conflnement to flve swimmer sizes, 2H = 5, and the 
periodicity at three times the conflnement, L = 3 x 2H. We will show later that this 
periodicity is large enough for the results to be independent of L because the walls 
prevent formation of flow structures larger than the conflnement. In the dilute theory 
described previously, it was assumed that the swimmers formed a layer at each wall. 
The formation of layers close to the conflning wall in the simulations is illustrated in 
flgure [5} which shows the swimmer concentration proflles as functions of X3 at different 
eflFective volume fractions. Results for both pushers and pullers are shown, as well as 
for swimmers where hydrodynamic effects are completely neglected. For very dilute 
systems, the results in all cases are fairly similar - layers form near each wall and there 
is a nearly zero concentration at small distances from the walls due to steric exclusion of 
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beads from the walls. However, the structure within the layer does depend on whether 
hydrodynamic interactions are included. With hydrodynamic interactions, no long- 
range orientational order is observed in the layer. Without hydrodynamic interactions, 
a two-dimensional nematic state is observed, similar to the three-dimensional nematic 
state observed in an unconfined system without hydrodynamic interaction ^12]. Once 
the concentration is high enough that the layer close to the wall is saturated, the 
swimmers form additional layers indicated by the secondary peaks in the profiles at 
high concentrations. 

Figure [6] shows the mean-squared-displacement (MSD) in the periodic (xi, X2) plane 
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as a function of time, for swimmers and fluid tracers, at two concentrations: (/)e = 0.1 
(^^e = 0.375) and (/)e = 0.3 {ipe = 1.125). At short times, the MSD is baUistic, reflecting 
the straight-hne motions of an isolated swimmer. The duration of this baUistic regime 
decreases as the concentration increases. At longer times the behavior becomes diffusive. 
Since Brownian motion is absent, the origin of this diflFusive regime is the interactions 
between the swimmers. Figure [7] shows the long-time diffusion coefficients Dg and as 
a function of the effective volume fraction and area fraction. At low concentrations the 
swimmers have a high effective diffusion coefficient reflecting the weak hydrodynamic 
interactions between them. The flow is only disturbed by a small number of swimmers 
so tracers diffuse very slowly. As the concentration is increased the diffusivity of the 
swimmers decreases, as the natural ballistic trajectories are increasingly perturbed by 
hydrodynamic interactions and collisions with other swimmers. Correspondingly, the 
naturally motionless tracers feel the motion of increasingly more swimmers and their 
diffusivity increases. 

According to the results in flgure [7[ the diffusion coefficient for both swimmers and 
tracers follow the dilute theory scalings { ^ ^ Dt ^ i/j^ ) developed in section [3] 
even at moderate concentrations. Note the difference in diffusivity and density proflles 
compared to our previous work [8j . For the volume fractions shown in the present work, 
we do not see the same dramatic increase in diffusivity and shift in the density proflles 
to the center of the channel. These differences seem to be due to the use of regularized 
forces in the present work, versus point forces in the former, and an excluded volume 
potential as shown in flgure 1. With this excluded volume potential and swimmer aspect 
ratio, obtaining much larger volume fractions is not possible because there is a largest 
volume fraction corresponding to the close-packed state. 

From section [3j we saw that the diffusivities in flgure [7] are related to the typical 
velocities of the swimmers and tracers. Figure [s] shows {v'^) and ('^g ) — 1 as functions of i/Jq. 
Both display an approximately linear increase with concentration, as predicted by the 
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Figure 7. Diffusion coefficients as a function of the effective volume fraction (top 
axis) and effective area fraction (bottom axis) for swimmers and tracers with L = 15 
and 2iJ 5. 



10 



A 



10" 





^ ^puller 








/, ^^___pusher 























10 



10 



10 



10 



Figure 8. Mean-squared velocity for swimmers, (y^^ — 1, and tracers, (yf) as a 
function of the effective area fraction for L = 15 and 2E. = 5. 



simple theory. We believe the deviation at small concentration is because of statistical 
errors due to the small number of swimmers in the domain at small concentrations. As 
with the diffusivities, the dilute scaling seems to hold even at concentrations well above 
ipQ = 0.1 for pushers. However, the puller simulations deviate from the dilute scaling at 
i/Jq > 0.2. The cause for this difference between pushers and pullers is unknown. 

The final property we examine for this initial confinement is the fiuid fiow generated 
by the swimmers. Figure [9] shows snapshots of the fiuid velocity field after more than 
1000 dimensionless times in a plane near the wall of the channel (X3 = 0.8H) for different 
concentrations cpe = 0.01,0.1 and 0.3 {ipe = 0.0384,0.375 and 1.125). These velocity 
fields are similar to those observed in some experiments ^J. Snapshots suggest that 
the size of the fiow structures is related to the wall separation 2H (as we will address 
later). We also examined the snapshots for three different locations in the channel, 
Xs = 0, 0.5i7, and 0.8i7, which are shown in figure 10 for (/)e = 0.1 {ipe = 0.375). 




Figure 9. Snapshots of the velocity field at xs = 0.8H {xs = H — 1/2) with L = 15 
and 2H = 5. (a) 0e = 0.01 (?Ae = 0.0383) (b) = 0.1 (V^e = 0.375) (c) 0e = 0.3 

(V^e = 1.125) 




Figure 10. Snapshots of the velocity field for = 0.1 {i/jq = 0.375) with L = 15 and 
2H = 5. (a) xs = 0.8H {xs = H - 1/2) (b) xs = 0.5H (c) xs = 
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Figure 11. Diffusion coefficients vs effective area fraction for swimmers and tracers 
for different confinements. For 2H = 5 and 2H = 10 we used L = 3 x 2H^ while for 
the monolayer 2H = 1 we used L = 20 x 2H. 



We now turn to the effect of the degree of confinement on the results. Simulations 
with a wider gap {2H = 10 and L = 3 x 2H) will be examined as well as a highly 
confined case of 2H = 1, L = 20 x 2i/, and with the swimmers restricted to lie in the 
midplane of the channel - a monolayer. Figures [TT] and [12] show the long-time diffusion 
coefficients and D^^ the mean-squared swimmer velocity deviation {v'^) — 1, and the 
mean-squared fluid velocity (v^) for all three confinements. The first observation is that 
all confinements follow the dilute theory scalings at low concentrations. In particular, 
we see that the swimmer diffusivities for 2H = 5 and 2H = 10 confinements collapse 
reasonably well when the concentration is represented as an effective area fraction. This 
supports the simple theory that assumes the swimmers form a layer next to each wall 
in the dilute limit. The swimmer diffusivities in the monolayer have the same scaling, 
though a different prefactor. This difference is most likely due to either a difference 
in the hydrodynamic interactions between swimmers in the monolayer because of the 
nearby walls or the inability of the swimmers to escape from the layer. 

The tracer diffusivities, while following the simple scaling with area fraction, do not 
collapse with ip^ with changing confinement. Because the tracers are not restricted to 
planes near the walls, at weak confinement (large H) the tracers can sample the lower 
velocities in the center of the channel, away from the swimmer layers near the walls. 
This difference results in a lower velocity of tracers with weaker confinement, as seen in 



figure 12 



From the earlier theoretical analysis, we found that uncorrelated swimmers produce 
swirls in the fluid, resulting in negative fluid correlations. The size of these negative 
regions were related to the separation of the walls, the size of a swimmer, and the 
separation of a swimmer from a wall. We found in the simulations swirls in the fluid as 
shown in figures [9] and 10, To make a quantitative comparison between the simulations 
and the theory of uncorrelated swimmers, we compare the fluid autocorrelation functions 
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Figure 12. Mean-squared velocity for swimmers, ('Ug ) — 1, and tracers, (^v^) vs effective 
area fraction. For 2H = 5 and 2H = 10 we used L = 3 x while for the monolayer 
2iJ = 1 we used L = 20 x 2H. 



in planes parallel to the slit. Figure [13] shows the velocity autocorrelation function at 
two X3-planes, X3 = 0.0 and 0.8i7, for 2H = 5 and 2H = 10, both with L = 3 x 2H and 
= 0.05. These results indicate, as foreshadowed by the velocity fields, that the scale of 
the velocity autocorrelation is set by the distance between the walls. The fluid velocities 
are correlated at short length scales and, before they become decorrelated, there is 
an anti-correlation. While these types of negative correlations might be attributed to 
swirls caused by collective behavior, it is necessary to show an explicit comparison 
with the uncorrelated swimmer response because independent swimmers can produce 



negative correlations due to the change in hydrodynamics caused by the wall. Figure 13 



shows such a comparison. First consider the comparison between curves in the plane 
Xs = 0.8 H. We see that the full simulation correlations are similar to the independent 
swimmer calculation. The difference between the curves in slightly larger in the plane 
X3 = 0. Recall that the independent swimmers are all placed at the walls, which produce 
fluid correlations in the center of the channel. The full simulations however also have 
some swimmers in the center of the channel. This additional contribution from the 
swimmers in the center may be sufficient to account for the difference. 

We now return to the issue of system size effects and our choice of L in the periodic 
directions. In order to explore the system size effects on the results, the confinement 
of 2H = 5 was kept constant while the periodicity of the box was increased from 
L = 3 X 2// to L = 4 X 2//, 5 X 2//, 6 X 2// and 12 X 2H. UnderhiU et al [12\ 
found that the tracer diffusivity in a three-dimensional periodic unconfined system has 
a dependence on system size. Because the tracer diffusivity in the unconfined system 
had a stronger system size dependence than the swimmer diffusivity and the velocities, 
we show in figure 14 the system size dependence of the tracer diffusion coefficient in the 
confined domain. No system size dependence is observed. We also did not observe any 
system size dependence of the swimmer diffusivity or swimmer and tracer velocities. 
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Figure 13. Fluid velocity autocorrelation function at (j)Q = 0.05 (a) 2H = 5 
(^e = 0.19) (b) 2H = 10 (^e = 0.37) 
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Figure 14. System size dependence of the tracer (fluid) diffusion coefficient at three 
effective volume fractions. 



This explicitly justifies our use of L = 3 x 2H in the data presented earlier. 

Finally, figure [15] illustrates the velocity autocorrelation at 2H = 5 for three system 
sizes L = 3 X 2i7, 6 x 2H and 12 x 2H. The curves have a slight change with system 
size at small separations, though they seem to converge quickly with increasing system 
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size. Clearly, the lack of changes with a system size in the confined system compared 
to the unconfined system is due to the introduction of an additional length scale once 
the system is confined. This is not surprising because the confinement restricts fiow 
structures larger than the confinement and modifies the hydrodynamic interactions 
between swimmers. If the collective behavior seen in the unconfined simulations is due to 
the long-range nature of the hydrodynamic interactions, screening of these interactions 
by the walls may affect the collective structures formed. The effect of hydrodynamic 
screening on the dynamics of polymers has been studied in detail [20l IHU [32] . In the 
polymer literature, it is the rate of decay of the fiuid disturbance with distance (as 
well as the rotational symmetry) that determines if screening is present or not. At 
large enough distance from a swimmer in the confined domain, the velocity disturbance 
decays as r~^. It is not clear if this decay is responsible for the results seen here - that 
the fiuid structures do not depend on system size provided the periodicity is larger than 
the scale of confinement. 
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5. Conclusions 

Simulations and theoretical analysis were used to study the dynamics of confined 
suspensions of self-propelled particles. The swimmers interact with one another via 
an excluded volume potential and hydrodynamic interactions through the fluid. These 
hydrodynamic interactions are altered by the conflning walls. 

We developed a simple theory of the motion of swimmers and tracers which captures 
the scalings of the diffusivity and velocities from the simulations. The theory assumes 
that in the dilute limit the swimmers form layers near the walls and execute a two- 
dimensional random walk within the layer. The theory also shows that even independent 
swimmers, with no collective behavior, can produce spatial fluid correlations in a 
conflned domain that do not appear in an unconflned domain. This is particularly 
important because swirls in the fluid and negative correlations have been cited in 
previous studies as evidence of collective behavior. Using this theory, we were able to 
better understand how the scale of these negative correlations in the simulations change 
with the degree of conflnement. In particular, the negative correlations in the center of 
the slit are governed by the separation of the walls, while the negative correlations near 
the walls are governed by the size of a swimmer and the separation of the swimmer from 
the wall. This is important for understanding experimental observations of swimming 
suspensions because many experiments are performed in the presence of walls. 
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